function hfile = smp_Smooth(hfile, srf, niter, mapno)
% SMP::Smooth  - smooth SMP map with neighbors of an SRF
%
% FORMAT:       [smp = ] smp.Smooth(srf [, niter [, mapno]]);
%
% Input fields:
%
%       srf         SRF object with required neighbor information
%       niter       number of iterations (default: 1)
%       mapno       which map number (default: 1)
%
% Output fields:
%
%       smp         object with one added, smoothed map

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'smp') || ...
    numel(srf) ~= 1 || ...
   ~isBVQXfile(srf, 'srf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
sbc = bvqxfile_getcont(srf.L);
if bc.NrOfVertices ~= sbc.NrOfVertices
    error( ...
        'BVQXfile:BadArgument', ...
        'SMP and SRF objects must match in NrOfVertices property.' ...
    );
end
if isempty(bc.Map)
    error( ...
        'BVQXfile:BadArgument', ...
        'SMP does not contain any maps.' ...
    );
end
if nargin < 3 || ...
   ~isa(niter, 'double') || ...
    numel(niter) ~= 1 || ...
    isinf(niter) || ...
    isnan(niter) || ...
    niter < 1 || ...
    niter > 10000
    niter = 1;
else
    niter = round(niter);
end
if nargin < 4 || ...
   ~isa(mapno, 'double') || ...
    numel(mapno) ~= 1 || ...
    isinf(mapno) || ...
    isnan(mapno) || ...
    mapno < 1 || ...
    mapno > numel(bc.Map)
    mapno = 1;
else
    mapno = round(mapno);
end

% get map values and neighbors
map = double(bc.Map(mapno).SMPData);
nei = sbc.Neighbors(:,2);

% prepare neighbors list
for vc = 1:numel(map)
    nei{vc}(end+1) = vc;
end

% iterate
for sc = 1:niter
    
    % create empty buffer for smoothed data
    mpb = zsz(map);
    
    % iterate over neighbors list
    for vc = 1:numel(mpb)
        mpb(nei{vc}) = mpb(nei{vc}) + map(vc) ./ numel(nei{vc});
    end
    
    % replace map
    map = mpb;
end

% add to object
bc.Map(end+1) = bc.Map(mapno);
if isempty(regexp(bc.Map(end).Name, ' sm\d+$'))
    bc.Map(end).Name = [bc.Map(end).Name sprintf(' sm%d', niter)];
else
    bc.Map(end).Name = regexprep(bc.Map(end).Name, ' sm\d+$', ...
        sprintf(' sm%d', niter));
end
bc.Map(end).SMPData = single(map);
bc.NrOfMaps = numel(bc.Map);

% set back in global storage array
bvqxfile_setcont(hfile.L, bc);
